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Abstract 



Biophysical models are often populated by a large number of input parameters that are difficult to 
predict or measure experimentally. The validity and robustness of a given model can be evaluated by a 
sensitivity test to its input parameters. In this study, we performed local (based on a Taylor-like method) 
Q^j and global sensitivity (based on Monte Carlo filtering techniques) analyses of a previously derived PAR1- 

mcdiatcd activation model of endothelial cells. This activation model previously demonstrated that 



peptide-activated PARI has a different receptor/G-protein binding affinity that favors Ga q activation over 
G Q! i2/i3 by approximately 800-fold. Interestingly, the present study shows that the parameter regulating 
the binding rate of activated PARI to Ga 12 /i3 is indeed important to obtain the expected RhoGTP 
response. Moreover, we show that the parameters representing the rate of PARI deactivation and the 
rate of PARI binding to G q , are the most important parameters in the system. Finally, we illustrate 
that the kinetic model considered in this study is robust and we provide complementary insights into the 
' biological meaning and importance of its kinetic parameters. 

Introduction 

Mathematical models of biophysical phenomena often involve a large number of input parameters [IHS] , 
such as reaction rates and initial concentrations, that must be either measured experimentally or inferred 
from similar biological systems. Moreover, experimentally measured parameters carry uncertainties due to 
experimental limitations, statistical analysis and different experimental conditions. One of the major goals 
in systems biology is to estimate how sensitive a computational model is to variations of its parameters 
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and study the effect of the parameter uncertainties on the model response. Sensitivity analysis aims 
at determining which parameters have the most influence on a predicted system behavior [BH5]. When 
the notion of "influence" is made quantitatively precise, sensitivity analysis can constitute a reliable 
robustness test for computational models [TO] . 

In this study, we performed sensitivity analysis on the mathematical model of PARl-mediated activa- 
tion of endothelial cells of [4] . A schematic representation of the main signaling pathways analyzed in [4] 
is shown in Fig.Q] A specific PARI agonist, SFLLRN, simultaneously activates two classes of G proteins: 
G q and G12/13 [11]. The a subunit of the G q protein activates the /3 isoforms of phospholipase C (PLC/3), 
which hydrolyze the lipid phosphatidylinositol 4,5-bisphosphate (PIP2) to generate the second messengers 
inositol 1,4,5-trisphosphatc (IP3) and diacylglycerol (DAG) |12j . The second messenger IP3 regulates the 
activity of the inositol trisphosphate receptors (IP3R) on the surface of the endoplasmic reticulum (ER), 
allowing the rapid release of Ca 2+ into the cytoplasm [T3"lll4j . Simultaneously, the a subunit of the G12/13 
protein activates the small GTPase RhoA, known to promote cytoskeletal changes [T5j . 

The signaling pathways described above are also common to a number of transduction systems me- 
diated by G-protein coupled receptors in other cell types. We have recently published a model of PARI 
signaling in platelets ( [13.:). 

In [3], several kinetic parameters had to be inferred because they were unknown and/or difficult to 
measure. In this work, we analyzed the influence of the chosen parameters to the system's output. We 
performed a sensitivity test based on the Taylor expansion of the output functions around the chosen 
input parameters (local analysis) and a sensitivity test based on Monte Carlo sampling techniques of the 
system parameters (global analysis) . Our results show that although the two techniques implemented are 
conceptually different, they lead, at least qualitatively, to consistent results. We show that the system 
is not sensitive to the majority of the parameters chosen in [3] and we identify three important nodes in 
the model. These results represent a further test for the validity of the model in [3] . Additionally, this 
analysis will help design more refined computational models, and will have the ultimate goal of identifying 
influential parameters and key signaling nodes with possible applications in biology and medicine. 
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Figure 1. Pictorial representation of PARl-mediated signaling pathways in endothelial 
cells. Upon cell stimulation by a PARI agonist peptide such as SFLLRN, PARI simultaneously 
activates the two classes of G proteins, G q and G12/13. The a subunit of the G q protein activates the j3 
isoforms of phospholipase C (PLC/3), which hydrolyze the lipid phosphatidylinositol 4,5-bisphosphate 
(PIP2) to generate the second messengers inositol 1,4,5-trisphosphate (IP3). IP3 then regulates the 
activity of the inositol trisphosphate receptors (IP3R) on the surface of the endoplasmic reticulum (ER) 
allowing the rapid release of Ca 2+ into the cytoplasm. Simultaneously, the a subunit of the G12/13 
protein activates the small GTPase Rho to promote cytoskeletal changes. Reproduced from g]. 

Methods 

Sensitivity Analysis by Taylor Expansion 

Consider a mathematical model defined by the following g-dimensional system of ordinary differential 
equations (ODEs) with initial conditions y D , 

[y -/<**) (1) 
y. =y(o), 



where y = (j/i, . . . , y q ) is the g-dimcnsional vector of system states and k = (ki, . . . , k n ) is the n- 
dimcnsional vector of input parameters. A simple estimate of the effect of the uncertainties in k on y, 
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Table 1. The sensitivity coefficients calculated by the Taylor-like formula at the peak time 
t; for Ca 2+ and RhoGTP 







RhoGTP 




k 


s+ 


s_ 


s+ 


S_ 


FUNCTIONAL ROLE 


1 


0.000 


0.002 


0.012 


0.038 


PARI activation by SFLLRN 


2 


0.000 


0.000 


0.000 


0.000 




3 


-0.001 


0.001 


-0.004 


-0.018 




7 


-0.006 


-0.006 


-0.193 


-0.197 


PARI deactivation 


8 


0.006 


0.006 


-0.051 


-0.050 


PARI binding Gq 


9 


-0.001 


-0.001 


0.008 


0.009 




10 


0.008 


0.010 


0.003 


-0.007 


Gq releasing GDP 


11 


0.000 


0.000 


0.000 


0.000 




12 


0.000 


0.000 


0.001 


0.001 


Gq binding GTP 


13 


0.000 


0.000 


0.000 


0.000 




14 


0.019 


0.024 


0.056 


0.053 


Gq releasing (3y 


15 


0.000 


0.000 


0.000 


0.000 




16 


0.000 


0.000 


0.000 


0.000 


Gq deactivation 


17 


0.000 


0.000 


0.014 


0.000 


PLC/3 binding GqGTP 


18 


0.000 


0.000 


0.000 


0.000 




19 


-0.040 


-0.034 


-0.018 


-0.003 


PLC/3 hydrolizing GqGTP 


20 


0.000 


0.000 


0.000 


0.000 


PLC/3 releasing GqGDP 


21 


0.000 


0.000 


0.000 


0.000 




22 


0.032 


0.042 


0.007 


0.013 


PLC/3 binding PIP2 


23 


0.000 


0.000 


0.000 


0.000 




24 


0.000 


0.000 


0.000 


0.000 


PLC/3 hydrolizing PIP2 


25 


-0.065 


-0.066 


-0.016 


-0.009 


Consumption of IP3 


26 


0.000 


0.000 


0.245 


0.247 


PARI binding G13 


27 


0.000 


0.000 


-0.041 


-0.041 




28 


0.000 


0.000 


0.039 


0.042 


G13 releasing GDP 


29 


0.000 


0.000 


0.000 


0.000 




30 


0.000 


0.000 


0.000 


0.000 


G13 binding GTP 


31 


0.000 


0.000 


0.000 


0.000 




32 


0.000 


0.000 


-0.001 


-0.001 


G13 releasing ^7 


33 


0.000 


0.000 


0.000 


0.000 




34 


0.000 


0.000 


-0.022 


-0.022 


G13 deactivation 


35 


0.000 


0.000 


0.106 


0.112 


GEF binding G13GTP 


36 


0.000 


0.000 


-0.107 


-0.111 




37 


0.000 


0.000 


-0.054 


-0.054 


GEF hydrolizing G13GTP 


38 


0.000 


0.000 


0.000 


0.000 


GEF releasing G13GDP 


39 


0.000 


0.000 


0.000 


0.000 




40 


0.000 


0.000 


0.249 


0.251 


GEFG13GTP binding RhoGDP 


41 


0.000 


0.000 


-0.017 


-0.017 




42 


0.000 


0.000 


0.017 


0.018 


Rho activation 


43 


0.000 


0.000 


-0.199 


-0.203 


R.ho deactivation 



The sensitiviy coefficients S%j{t, k) 



d\n( yi ) 
81n(fc 3 -) 



calculated at the nominal peak time t; for = Ca 2+ 



k=k 



and variations Ak = ±5% (columns two and three) and yi = RhoGTP and variations Ak = ±5% 
(columns four and five). Sensitivity coefficients larger than 10% are in bold font. The nominal peak 
time ti is defined in ([5]) and the coefficients SV, are defined in ([3]). 



that is, the variation of y caused by a variation Ak in k is 



Ay(t,k)=y(t,k + Ak)-y(t,k), 



where y(t, k) is a solution of (JTJ at time t with parameters k = k. Assuming Ak <C 1, by the Taylor's 
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expansion in the variables k about k, 



Ay(t, k) = D(t, k)Ak + 0(| Ak| 2 ) 



(2) 



where 0(|Ak| 2 ) is a (/-dimensional vector infinitesimal of order no less than 2 with respect |Ak|, and 
D(t, k) is the n x q matrix of the partial derivatives at k of entries 



i = l,...,q 
k=k j = l,...,n 



Then one takes a dimcnsionlcss version of D(i, k) as a measure of the sensitivity of the system at time t, 
with respect to the set of parameters k. Precisely one introduces a sensitivity matrix S(t, k) of entries 



S«(t,k) 



dln(yi 



<91n(%) 



k 3 dyi(t,k) 



k=k 



j/i(t,k) 



(3) 



k=k 



and takes Sij(t, k) as a measure of the sensitivity of the system state yi(t, k) at time t with respect to the 
kinetic parameter kj about the nominal values k. In the context of the kinetic model in [4] we selected 
the states 

2/ 4 (<,k) = [Ca 2 +](i,k) 

(4) 

2/i(i,k) = [RhoGTP](i,k) 

and computed SV,(ti,k) at time t; where the nominal [Ca 2+ ](i, k) and [RhoGTP](i, k) attain their max- 
imum values, e.g., 

2/i(t;,k) = max?/i(t,k). (5) 

The sensitivity coefficients 5ij(ti,k) for Ca 2+ are reported in Table Q] in columns two and three along 
those for RhoGTP in columns four and five. The partial derivatives dyi/dkj were approximated by 
discrete differences with values of ±5% of the nominal values and indicated respectively with the symbols 
S + and S-. 

The total output of y(-,k) over the average time course T of the experiment, is 



z(k) 



y(t,k)dt. 



(6) 



G 



Table 2. The sensitivity coefficients calculated by the Taylor-like formula over the course 
of the entire simulation T = 600 s for Ca 2+ and RhoGTP 







RhoGTP 




k 


£+ 


£_ 


£+ 


s_ 


FUNCTIONAL ROLE 


1 


0.000 


0.000 


0.008 


-0.017 


PARI activation by SFLLRN 


2 


0.000 


0.000 


0.000 


0.000 




3 


0.000 


0.000 


-0.003 


-0.013 




7 


0.000 


0.000 


-0.449 


-0.443 


PARI deactivation 


8 


0.000 


0.000 


0.021 


0.040 


PARI binding Gq 


9 


0.000 


0.000 


-0.007 


-0.007 




10 


0.000 


0.000 


-0.001 


-0.012 


Gq releasing GDP 


11 


0.000 


0.000 


0.000 


0.000 




12 


0.000 


0.000 


-0.001 


-0.001 


Gq binding GTP 


13 


0.000 


0.000 


0.000 


0.000 




14 


0.000 


0.000 


-0.006 


-0.016 


Gq releasing (3y 


15 


0.000 


0.000 


0.000 


0.000 




16 


0.000 


0.000 


0.000 


0.000 


Gq deactivation 


17 


0.000 


0.000 


0.011 


0.000 


PLC/3 binding GqGTP 


18 


0.000 


0.000 


0.000 


0.000 




19 


0.000 


0.000 


-0.022 


0.002 


PLC/3 hydrolizing GqGTP 


20 


0.000 


0.000 


0.000 


0.000 


PLC/3 releasing GqGDP 


21 


0.000 


0.000 


0.000 


0.000 




22 


0.000 


0.000 


0.001 


0.019 


PLC/3 binding PIP2 


23 


0.000 


0.000 


0.000 


0.000 




24 


0.000 


0.000 


0.000 


0.000 


PLC/3 hydrolizing PIP2 


25 


-0.051 


-0.057 


-0.009 


0.031 


Consumption of IP3 


26 


0.000 


0.000 


0.423 


0.429 


PARI binding G13 


27 


0.000 


0.000 


-0.070 


-0.071 




28 


0.000 


0.000 


0.068 


0.074 


G13 releasing GDP 


29 


0.000 


0.000 


0.000 


0.000 




30 


0.000 


0.000 


0.000 


0.000 


G13 binding GTP 


31 


0.000 


0.000 


0.000 


0.000 




32 


0.000 


0.000 


0.000 


0.000 


G13 releasing ^7 


33 


0.000 


0.000 


0.000 


0.000 




34 


0.000 


0.000 


-0.199 


-0.202 


G13 deactivation 


35 


0.000 


0.000 


0.064 


0.069 


GEF binding G13GTP 


36 


0.000 


0.000 


-0.066 


-0.066 




37 


0.000 


0.000 


-0.488 


-0.504 


GEF hydrolizing G13GTP 


38 


0.000 


0.000 


0.000 


0.000 


GEF releasing G13GDP 


39 


0.000 


0.000 


0.000 


0.000 




40 


0.000 


0.000 


0.429 


0.434 


GEFG13GTP binding RhoGDP 


41 


0.000 


0.000 


-0.030 


-0.030 




42 


0.000 


0.000 


0.029 


0.031 


Rho activation 


43 


0.000 


0.000 


-0.402 


-0.421 


Rho deactivation 



The sensitiviy coefficients S,j (k) 



dz 



Zi(k) dk 3 



with z(k) = y(i, k)dt for j/j = Ca 2+ and 



k=k 



variations Ak = ±5% (columns two and three) and yi = RhoGTP and variations Ak = ±5% (columns 
four and five). Sensitivity coefficients larger than 10% are in bold font. 



For the states in (0| we chose in our simulations T = 600s. The vector z(k) is independent of time and it 
can be expanded in Taylor's series with respect to k, as in @ with D(i, k) replaced by its time integral 
over (0,T). The sensitivity coefficients of z(k) about the nominal values k are 



Ey(k) 



(7) 



k=k 
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The sensitivity coefficients for Ca 2+ are reported in Table [2] in columns two and three along those 
for RhoGTP in columns four and five. The partial derivatives dyi/dkj were approximated by discrete 
differences with values of ±5% of the nominal values and indicated respectively with the symbols E + and 
E_. 

The nominal values of the parameters kj used in the sensitivity analysis described above are given in 
Table 1 

Monte Carlo Filtering Sensitivity Analysis 

The method is based on estimating the uncertainties distributions p(k) = {pi(k\), . . . ,p n (k n )} of the 
parameters k = (k\, . . . , k n ), each ranging over the intervals 

k 3 , 6 (kf\ fcf ax ), j = l,...,n (8) 

and cumulatively generating the probability measure 

p(k)dk = pi(ki) ■ ■ -p„(fc„) dki ■ ■ ■ dk n 

over the space of parameters 

n 

n (kf n ,kf ax ). 

If these ranges and uncertainties distributions were known, regarding each of the yi(t,h) as a random 
variable depending on the random k one computes the mean 



and, the variance 



where 



(yM) = J w (t,k) P (k)dk 



of(t) = (yUt)) ~ (Viit)) 2 , 



(y?(*)> = /•••/! w?(t,k) P (k)dk. 



Table 3. Symbolic reaction schemes and effective kinetic parameters used in the simulations of the PARl-mediated 
activation model of endothelial cells proposed in [4]. 



Symbolic Reactions 


k_> | k<_ | References 


Reactions govern-in 


} PARI activation 


Agonist + PARI ^ PARI* 


6.00E+04 M _1 s _1 {ki) 


1.00E-03 s _1 (fc 2 ) 


16 


Agonist — > null 


2.00E-01 s _1 (fc 3 ) 




□ 


PARI* -> null 


2.00E+01 s _1 (fe 7 ) 




□ 


Reactions governing Gg activation 


PARI* + GqGDP ■ 07 ^ PARI* • G q GDP ■ 0i 


1.00E+08 M _1 s _1 (fe 8 ) 


1.00 s _1 (fc 9 ) 


eh 


PARI* • GqGDP ■ 0i ^ PARI* ■ G q ■ 01 + GDP 


5.00 s _1 (fcm) 


1.00E+06 M _1 s _1 (fcn) 


eh 


PARI* • G q ■ 01 + GTP ^ PARI* ■ G q GTP • /3 7 


1.00E+06 M _1 s _I (fci 2 ) 


1.00E-01 s _1 (fei 3 ) 


eo 


PARI* • G q GTP ■ 01 ^ PARI* + G q GTP + 0i 


2.00 s _1 (fc i4 ) 


1.00E+07 M^s -1 (fci 5 ) 


E2 


GqGTP -> GqGDP 


2.00E-02 s _1 (fcie) 




fIM20l 


Reactions governing IPS and DAG generation 


PLC/3 + G q GTP ^ PLC0 ■ G q GTP 


5.00E+08 M _i s _i (fci 7 ) 


5.00 s" 1 (fe 18 ) 




PLC0 • GqGTP -> PLC/3 • GqGDP 


1.5E+01 s" 1 (fci 9 ) 




mm 


PLC,3 • GqGDP ^ PLC0 + GqGDP 


1.00E+05 s" 1 (fc 20 ) 


1.00E+02 M _i s _I (fc 2 i) 


22 


PLC/3 • GqGTP + PIP2 ^ PLC/3 ■ G q GTP • PIP2 


1.00E+09 M _i s _i (fe 22 ) 


1.00 s" 1 (fc 23 ) 


[22] 


PLC0 • GqGTP • PIP2 -> PLC0 • G q GTP + IP3/DAG 


1.00E+02 s" 1 (fc 24 ) 




22 


IP3 -> null 


2.40E-02 s" 1 (fc 25 ) 




® 


Reactions governing Ca2+ mobilization 


As described previously 1 1 4=1 


Reactions governing G12/13 activation 


PARI* + G 12/13 GDP ■ l3~f ^ PARI* ■ G 12/13 GDP ■ 0j 


1.00E+08 M _1 s _1 (fc 26 ) 


1.00 s _1 (k 27 ) 


|l7] 


PARI* • Gi 2/13 GDP ■ 0i ^ PARI* ■ G 12/13 ■ 01 + GDP 


6.00 s" 1 (fc 28 ) 


1.00E+06 M _1 s _1 (fc 29 ) 


ED 


PARI* ■ G 12/13 ■ 01 + GTP ^ PARI* ■ G 12 /i 3 GTP ■ 0-y 


1.00E+06 M _i s _1 (fc 30 ) 


1.00E-01 s _1 (fc 3 i) 


ED 


PARI* • G 12/13 GTP ■ 0i ^ PARI* + G 12/13 GTP + 0-y 


2.00 s- 1 (fe 32 ) 


1.00E+07 M-^s -1 (fc 33 ) 


E3 


Gi 2/ i 3 GTP -> G 12/13 GDP 


4.00E-03 s" 1 (fc 34 ) 




23 24 


Reactions governing Rho activation 


GEF Rho + G 12/13 GTP ^ GEFnho ■ G 12/13 GTP 


1.00E+09 M _1 s _1 (fc 35 ) 


3.00 s _1 (fc 36 ) 


[20] 


GEF Rho ■ G 12/13 GTP -> GEFRho ■ G 12/13 GDP 


1.00E-2 s _1 (fc 36 ) 




23 


GEF Rho ■ G 12/13 GDP ^ GEF Rho + G 12/13 GDP 


1.00E+06 s" 1 (fc 38 ) 


3.00 M _1 s _1 (/£ 39 ) 




GEF Rho ■ Gi 2/13 GTP + RhoGDP ^ GEF Rho • Gi 2/13 GTP ■ RhoGDP 


7.70E+09 M _1 s _1 (fc 4 o) 


7.70E+02 s" 1 (fc 4i ) 


s 


GEF Rho ■ Gi 2/13 GTP ■ RhoGDP -> GEF Rho ■ Gi 2 /i 3 GTP + RhoGTP 


1.03E+04 (fe 42 ) 






RhoGTP -> RhoGDP 


6.00E-02 s _1 (fc 43 ) 




[25] 



The nominal values of the parameters adopted in the PARl-mediated activation model of endothelial cells proposed in [4]. The kinetic 
constants are numbered sequentially starting from the first fc_j. (fci) and proceeding in reading order with the exception of reaction 
"PARI* — > null" which is regulated by the constant kr- 
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The variance is a measure of the influence, and relative importance, of the input k to the output y(t, k) 

(£21). 

In practice the process is implemented in a less quantitative way, by trading the information com- 
ing from the variances (Ji(t) with those coming from a biologically motivated objective functional, (fobj 
depending on one or several ?/i(i,k) |26j . Having chosen an objective functional g bj, one identifies a 
biologically acceptable range for the objective function g bj- For example for a threshold value ^ thres of 
the objective functional <? bj one might define [26] 

<7obj < g thlcs as the acceptable range 
ffobj > g thros as the non acceptable range. 

Then for each kj one determines the probability distributions f\ (kj ) and f 2 (kj ) of those values of kj that 
output the system in the acceptable or non acceptable range respectively. For each of these, and for each 
kj , calculate the cumulative frequency distributions 

cfi(kj) = — ^ tjtt, £ = 1,2 (10) 

max reiki) 

and calculate the Kolmogorov-Smirnov coefficients di t2 (kj) by the maximum distance function 

di >2 (k j )= sup \chikj)- cf 2 {kj)\. (11) 

These coefficients are taken as a measure of the relative importance of each parameter kj on the model 
output. The larger the value of di^ikj), the more important is kj in producing the pre-defined system 
output [6HS[26]. 

In the context of the kinetic model of [4] the method is implemented as follows: 

1. Select nominal values k = (ki, ■ . ■ , fc n ) as those originally adopted in the model of [1] (see Table [3]). 
As ranges in ([SJ we take intervals spanning from 1/10 to 10 times these nominal values, e.g., 

kj ;= lOAjj, and kj z:zz r^kj. 
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On a log 10 scale these are symmetric intervals about log 10 kj. 

2. Generate M random n-tuples of numbers 

k m = (fci.m, ■ • ■ , k n ,m) for m = 1, . . . , M 

by uniformly sampling the In kj from their respective symmetric uncertainty ranges defined above. In 
the simulations we used M = 10, 000. Different sampling distributions (e.g. gaussian, exponential) were 
seen not to qualitatively affect the results. 

3. Solve the system (JTJ) for each random choice of k rn , and compute the functions y(t, k m ). Solve also 
(JTJ for the nominal values k = k to get the functions y(i, k). Then for each % = 1, . . . , q introduce two 
kinds of objective output functions of the mth trial as follows: 

Objective Function at Times U: 

gobj ;J (ti;m) = [j/iOi.k™) - t/i(ti,k)] 2 (12) 

where m 6 {1, . . . , M} is the mi/i random trial described in the previous step and the times ti are de- 
fined in ©. For the states in (0}, this function measures the variations of the of [Ca 2+ ](i Ca 2+, k) and 
[RhoGTP](£p{ho7 k), from their nominal maximum values, (see Table 2] columns two and three respec- 
tively) . 

Objective Function for Total Time: 

Let z(k) be the total state output as in ^ and for i = 1, . . . , q set 

G obj;i (m) = [^(k m )-^(k)] 2 (13) 

where m £ {1, . . . , M} is the mth random trial. For the states in ((4]), the function G bj;i(w) measures 
the perturbation of the total outputs of Ca 2+ and RhoGTP over the time course T of the experiment, 
from their nominal values. The results are compared in Table [5] in columns two and three. 
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Table 4. The sensitivity coefficients calculated by Monte Carlo filtering at the nominal 
peak time t ; for Ca 2+ and RhoGTP 





di,2 




k 


Ca 2 + 


RhoGTP 


FUNCTIONAL ROLE 


1 


0.095 


0.084 


PARI activation by SFLLR.N 


2 


0.031 


0.026 




3 


0.113 


0.093 




7 


0.160 


0.285 


PARI deactivation 


8 


0.091 


0.153 


PARI binding Gq 


9 


0.037 


0.043 




10 


0.072 


0.017 


Gq releasing GDP 


11 


0.017 


0.022 




12 


0.018 


0.021 


Gq binding GTP 


13 


0.041 


0.019 




14 


0.153 


0.061 


Gq releasing {3y 


15 


0.023 


0.012 




16 


0.026 


0.024 


Gq deactivation 


17 


0.024 


0.033 


PLC/3 binding GqGTP 


18 


0.030 


0.038 




19 


0.421 


0.034 


PLC,3 hydrolizing GqGTP 


20 


0.042 


0.013 


PLC/3 releasing GqGDP 


21 


0.029 


0.029 




22 


0.394 


0.017 


PLC/3 binding PIP2 


23 


0.031 


0.030 




24 


0.034 


0.021 


PLC/3 hydrolizing PIP2 


25 


0.413 


0.018 


Consumption of IP3 


26 


0.035 


0.391 


PARI binding G13 


27 


0.016 


0.105 




28 


0.027 


0.068 


G13 releasing GDP 


29 


0.018 


0.026 




30 


0.042 


0.027 


G13 binding GTP 


31 


0.017 


0.015 




32 


0.038 


0.019 


G13 releasing /37 


33 


0.023 


0.048 




34 


0.017 


0.073 


G13 deactivation 


35 


0.021 


0.158 


GEF binding G13GTP 


36 


0.024 


0.196 




37 


0.026 


0.170 


GEF hydrolizing G13GTP 


38 


0.016 


0.043 


GEF releasing G13GDP 


39 


0.026 


0.025 




40 


0.042 


0.385 


GEFG13GTP binding RhoGDP 


41 


0.024 


0.104 




42 


0.029 


0.031 


Rho activation 


43 


0.023 


0.350 


Rho deactivation 



The sensitivity coefficients measured as Kolmogorov-Smirnov distances d\ t 2 estimated in the ranges 
[1/10, 10]fcj at the nominal peak time ti for yi = Ca 2+ (column two) and yi = RhoGTP (column three). 
Sensitivity coefficients larger than 10% are in bold font. The nominal peak time ti is defined in (0 and 
the objective function used in the simulations is defined in (|12p . 



4. Introduce threshold values 

1 M 

Mm=1 (14) 
M m=1 

The mth random trial and its parameters k m are deemed acceptable according to the criterion in §§§ for 
each of these objective functions and their respective threshold values. On the basis of this classification, 
generate the probability distribution functions fi(kj) and fiikj) of acceptable and unacceptable values, 
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Table 5. The sensitivity coefficients calculated by Monte Carlo filtering over the course of 
the entire simulation T = 600 s for Ca 2+ and RhoGTP 





di,2 




k 


Ca 2 + 


RhoGTP 


FUNCTIONAL ROLE 


1 


0.042 


0.032 


PARI activation by SFLLR.N 


2 


0.034 


0.029 




3 


0.045 


0.071 




7 


0.171 


0.208 


PARI deactivation 


8 


0.109 


0.049 


PARI binding Gq 


9 


0.041 


0.017 




10 


0.044 


0.017 


Gq releasing GDP 


11 


0.017 


0.018 




12 


0.046 


0.021 


Gq binding GTP 


13 


0.020 


0.029 




14 


0.034 


0.023 


Gq releasing {3y 


15 


0.030 


0.019 




16 


0.031 


0.038 


Gq deactivation 


17 


0.021 


0.050 


PLC/3 binding GqGTP 


18 


0.020 


0.028 




19 


0.157 


0.022 


PLC,3 hydrolizing GqGTP 


20 


0.035 


0.030 


PLC/3 releasing GqGDP 


21 


0.042 


0.011 




22 


0.085 


0.020 


PLC/3 binding PIP2 


23 


0.020 


0.013 




24 


0.028 


0.018 


PLC/3 hydrolizing PIP2 


25 


0.216 


0.043 


Consumption of IP3 


26 


0.032 


0.141 


PARI binding G13 


27 


0.034 


0.041 




28 


0.026 


0.044 


G13 releasing GDP 


29 


0.031 


0.048 




30 


0.033 


0.028 


G13 binding GTP 


31 


0.022 


0.024 




32 


0.050 


0.021 


G13 releasing /37 


33 


0.029 


0.017 




34 


0.038 


0.066 


G13 deactivation 


35 


0.034 


0.058 


GEF binding G13GTP 


36 


0.029 


0.111 




37 


0.033 


0.125 


GEF hydrolizing G13GTP 


38 


0.058 


0.011 


GEF releasing G13GDP 


39 


0.033 


0.016 




40 


0.044 


0.105 


GEFG13GTP binding RhoGDP 


41 


0.045 


0.052 




42 


0.052 


0.039 


Rho activation 


43 


0.037 


0.238 


R.ho deactivation 



The sensitivity coefficients measured as Kolmogorov-Smirnov distances d\ t % estimated in the ranges 
[1/10, 10] kj over the course of the entire simulation T = 600 s for yi = Ca 2+ (column two) and 
Hi = RhoGTP (column three). Sensitivity coefficients larger than 10% are in bold font. The objective 
function used in the simulations is defined in (|13l) . 



relative to g bj;i(ti; m ) an d G bj ; i(?7i) respectively. 

5. Calculate the cumulative frequency distributions cfe(kj) as in (fTU)) relative to each objective function, 
and the corresponding Kolmogorov-Smirnov coefficient d\_2{kj) as in (jlip . The larger the value of d\^{kj) 
the higher the sensitivity of the system to the variation of the corresponding parameter. 



All the calculations were performed on a MATLAB (R2009b, The Mathworks, Natick, MA) platform. 
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Times to Peak and Peak Values 

Table 6. The sensitivity coefficients calculated by Monte Carlo filtering for the peak 
values of Ca 2 + and RhoGTP 





di,2 




k 


Ca 2 + 


RhoGTP 


FUNCTIONAL ROLE 


1 


0.091 


0.082 


PARI activation by SFLLRN 


2 


0.030 


0.027 




3 


0.102 


0.100 




7 


0.218 


0.280 


PARI deactivation 


8 


0.160 


0.120 


PARI binding Gq 


9 


0.063 


0.035 




10 


0.073 


0.014 


Gq releasing GDP 


11 


0.026 


0.012 




12 


0.028 


0.026 


Gq binding GTP 


13 


0.054 


0.016 




14 


0.138 


0.055 


Gq releasing /3^y 


15 


0.024 


0.017 




16 


0.027 


0.021 


Gq deactivation 


17 


0.020 


0.036 


PLC/3 binding GqGTP 


18 


0.038 


0.047 




19 


0.404 


0.030 


PLC/3 hydrolizing GqGTP 


20 


0.029 


0.016 


PLC/3 releasing GqGDP 


21 


0.032 


0.035 




22 


0.438 


0.019 


PLC/3 binding PIP2 


23 


0.038 


0.019 




24 


0.014 


0.017 


PLC/3 hydrolizing PIP2 


25 


0.386 


0.022 


Consumption of IP3 


26 


0.041 


0.375 


PARI binding G13 


27 


0.017 


0.112 




28 


0.028 


0.059 


G13 releasing GDP 


29 


0.031 


0.018 




30 


0.022 


0.025 


G13 binding GTP 


31 


0.019 


0.021 




32 


0.021 


0.019 


G13 releasing /?7 


33 


0.018 


0.038 




34 


0.038 


0.077 


G13 deactivation 


35 


0.023 


0.160 


GEF binding G13GTP 


36 


0.027 


0.206 




37 


0.022 


0.130 


GEF hydrolizing G13GTP 


38 


0.036 


0.038 


GEF releasing G13GDP 


39 


0.024 


0.022 




40 


0.037 


0.370 


GEFG13GTP binding RhoGDP 


41 


0.021 


0.096 




42 


0.024 


0.035 


Rho activation 


43 


0.020 


0.344 


Rho deactivation 



The sensitivity coefficients measured as Kolmogorov-Smirnov distances dx,2 estimated in the ranges 
[1/10, 10] kj for the peak values of y% = Ca 2+ (column two) and yi = RhoGTP (column three). 
Sensitivity coefficients larger than 10% are in bold font. The objective function used in the simulations 
is defined in (fT6|) and the peak values are defined in ((T7)) . 



According to {[5J the times are those at when the nominal states (-, k) achieve their peak value. 
The states ?/.;(•, k) however achieve their peak values at times t^ eak which, in general differ from ti and 
are in general functions of k. Set 



2/r x (k) = W (t? eak (k) J k) - max^,k). 



(15) 
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Table 7. The sensitivity coefficients calculated by Monte Carlo filtering for the peak times 

tf ak of Ca 2 + and RhoGTP 





di,2 




k 


Ca 2 + 


RhoGTP 


FUNCTIONAL ROLE 


1 


0.075 


0.043 


PARI activation by SFLLRN 


2 


0.032 


0.028 




3 


0.031 


0.034 




7 


0.201 


0.530 


PARI deactivation 


8 


0.134 


0.441 


PARI binding Gq 


9 


0.063 


0.161 




10 


0.081 


0.076 


Gq releasing GDP 


11 


0.071 


0.032 




12 


0.074 


0.032 


Gq binding GTP 


13 


0.033 


0.045 




14 


0.199 


0.338 


Gq releasing {3y 


15 


0.067 


0.045 




16 


0.076 


0.023 


Gq deactivation 


17 


0.052 


0.042 


PLC/3 binding GqGTP 


18 


0.033 


0.032 




19 


0.148 


0.039 


PLC,3 hydrolizing GqGTP 


20 


0.025 


0.040 


PLC/3 releasing GqGDP 


21 


0.035 


0.046 




22 


0.324 


0.039 


PLC/3 binding PIP2 


23 


0.051 


0.040 




24 


0.066 


0.040 


PLC/3 hydrolizing PIP2 


25 


0.318 


0.022 


Consumption of IP3 


26 


0.038 


0.083 


PARI binding G13 


27 


0.023 


0.039 




28 


0.032 


0.059 


G13 releasing GDP 


29 


0.040 


0.030 




30 


0.027 


0.046 


G13 binding GTP 


31 


0.052 


0.016 




32 


0.056 


0.027 


G13 releasing /37 


33 


0.047 


0.065 




34 


0.067 


0.314 


G13 deactivation 


35 


0.044 


0.063 


GEF binding G13GTP 


36 


0.064 


0.083 




37 


0.034 


0.448 


GEF hydrolizing G13GTP 


38 


0.029 


0.044 


GEF releasing G13GDP 


39 


0.030 


0.042 




40 


0.067 


0.045 


GEFG13GTP binding RhoGDP 


41 


0.046 


0.051 




42 


0.060 


0.034 


Rho activation 


43 


0.029 


0.059 


R.ho deactivation 



The sensitivity coefficients measured as Kolmogorov-Smirnov distances di_2 estimated in the ranges 
[1/10, 10] kj for the peak times t? eak of yi — Ca 2+ (column two) and yi = RhoGTP (column three). 
Sensitivity coefficients larger than 10% are in bold font. The objective function used in the simulations 
is defined in ([16|) and the peak times tf eafc are defined in (fT5|) . 



Dose responses for different concentrations of agonists are measured at t^ oa . In practice, for a set of 
parameters k, which in general is unknown, this is done by recording the experimental outputs j/" lax (k). 
For small variations of k about its nominal vector k, the sensitivity of ax (k) could be theoretically 
"measured" as in @, by the sensitivity matrix of entries 

h ( %(tp k (k),k) %(tf ak (k) ; k)9tr aJt (k) ] 

k=& yf iax (k)l dkj ddt dk, J k=lc 



-.peakp--^ 9 In 



opcaK/, 

« 1 ' 91n(fc,) 
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This formula however requires the form of the functions k — > t^ eak (k), which are in general not know. 
Alternatively, the analysis can be carried by the filtering method, by introducing two new objective 
functions: 

F^ ak (m)=[y^(k m )-y i (ti,k)] 2 

(16) 

TP cak (m)= [tr k (k m )-i,] 
where the vector k m is the output of the mth random trial. Threshold values and distribution functions 
can be defined and determined as above. For yi = [Ca 2+ ] and y\ = [RhoGTP], denote their peak times 
by t^,° a + and t^j^ respectively, and set 

[Ca 2 +f ax (k) = [Ca 2+ l(tri(k),k) 
[RhoTGP] max (k) = [RhoGTP] (tP° h a k (k),k). 

Then the first of (fT6|) measures the variations of the peak values of [Ca 2+ ] and [RhoGRP] from their peak 
nominal values, whereas the second of (|16[) measures the variation of their times to peak t^° ak and t^^ 
respectively, from their nominal values t Ce 2+ and ipjioGTP- The results are reported in Tables [6] and [7] 
respectively. 



Results and Discussion 

We performed sensitivity analysis of a previously derived PARl-mediated activation model of endothelial 
cells. We used two different techniques, one based on the Taylor expansion of the system's output around 
the nominal values of its input parameters kj and the second based on Monte Carlo sampling techniques 
of the model parameters. The analysis based on the Taylor expansions @ and leading to the sensitivity 
matrices iSy(^,k) in §3§ and Sy(k) in ([7]), imposes two restrictions. The first is that |Ak| <C 1 with 
the notion of "smallness" depending on a predefined notion of smallness of |Ay|. The second is that 
these matrices measure the relative variation of j/i and with respect to kj , by keeping all the remaining 
parameters fixed, thereby neglecting the cumulative effects of Ak. 

The analysis based on the Kolmogorov-Smirnov test, does not require the range of k to be small, 
however it does require that the uncertainty distributions p(k) be known. The method leaves open the 
choice of the objective function. This on the one hand affords the flexibility in tailoring the objective 
function to specific experimental processes, and on the other hand allows non quantitative elements in 
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2 10 100 200 2 10 100 200 



Figure 2. Cumulative distribution functions of the most influential kinetic parameter k-j 
for the output j/j = Ca 2+ . The red curves indicate the distributions relative to the "acceptable" sets 
of parameters, whereas the blue ones denote those corresponding to the "unacceptable" set of 
parameters (see Methods). Panel A: Monte Carlo filtering method results at the nominal peak time tj. 
Panel B: Monte Carlo filtering results for the entire time-course of the simulation, T = 600 s. Panel C: 
Monte Carlo filtering results for the peak values [Ca 2+ ] max of the output function. Panel D: Monte 
Carlo filtering results for the peak times tP eak of the output function. The maximum distances between 
the acceptable and unacceptable distributions, are reported in Tables 0H7] respectively and are 
defined in (fTTj) . The nominal peak time tj is defined in flSj) and the objective functions used are 
respectively (fT2|) . (fT5)l , and ([To]) defined in Methods. 

the analysis. 

The two methods being complementary, we performed a sensitivity analysis by using both of them, 
on the states in Q arising from the mathematical model of [3J. 

First we analyzed the sensitivity of Ca 2+ and RhoGTP at their nominal peak values ((5]) to variations of 
the parameters k. Using the Taylor expansion method, we computed the sensitivity coefficients $y(ti, k) 
introduced in ([3]) and reported them in Table [TJ Using the Monte Carlo filtering method, starting from 
the objective function g^^^t^m) introduced in (|12[) , we computed the Kolmogorov-Smirnov coefficients 
by the distance function in (fTTj) , and reported them in Table [4] 
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Figure 3. Cumulative distribution functions of the most influential kinetic parameter k-j 
for the output j/j = RhoGTP. The red curves indicate the distributions relative to the "acceptable" 
sets of parameters, whereas the blue ones denote those corresponding to the "unacceptable" set of 
parameters (see Methods). Panel A: Monte Carlo filtering method results at the nominal peak time tj. 
Panel B: Monte Carlo filtering results for the entire time-course of the simulation, T = 600 s. Panel C: 
Monte Carlo filtering results for the peak values [RhoGTP] max of the output function. Panel D: Monte 
Carlo filtering results for the peak times tP eak of the output function. The maximum distances between 
the acceptable and unacceptable distributions, di j2 are reported in Tables 0H7] respectively and are 
defined in (fTTj) . The nominal peak time tj is defined in ([S]) and the objective functions used are 
respectively (fT2|) . (fT3")) , and (JTHJ) defined in Methods. 

Then we analyzed how parameters fluctuations affect the total response of Ca 2+ and RhoGTP over 
the whole time course T of the experiments, quantified by the integrated states ([6]). We computed first 
the sensitivity coefficients Sy(k) in (J7J, and reported in Table [5J Then, starting from the objective 
function G bj;i(T-) introduced in (|13|). we computed the Kolmogorov-Smirnov coefficients as in (|lip and 
reported in Table [5j 

Finally, by the same Monte Carlo filtering procedure we investigated the sensitivity of the peak values 
of Ca 2+ and RhoGTP and their relative times to peak for cumulative variations of all the parameters kj in 
intervals of 2 orders of magnitude with respect to kj. The results are in Tableland Table [7] respectively. 
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le+07 le+08 le+09 le+07 le+08 le+09 



k 26 k 26 

Figure 4. Cumulative distribution functions of the influential kinetic parameter k-26 for the 
output yi = RhoGTP. The red curves indicate the distributions relative to the "acceptable" sets of 
parameters, whereas the blue ones denote those corresponding to the "unacceptable" set of parameters 
(see Methods). Panel A: Monte Carlo filtering method results at the nominal peak time t;. Panel B: 
Monte Carlo filtering results for the entire time-course of the simulation, T = 600 s. Panel C: Monte 
Carlo filtering results for the peak values [RhoGTP] max of the output function. Panel D: Monte Carlo 
filtering results for the peak times t? eak of the output function. The maximum distances between the 
acceptable and unacceptable distributions, dip are reported in Tables [4HZ1 respectively and are defined 
in (jlip . The nominal peak time t; is defined in ([5]) and the objective functions used are respectively 
([12]) . (fT3|) . and ^T6j) defined in Methods. 

As a way of analyzing and comparing the results in these tables, we deemed a parameter impor- 
tant for a given state ?/i(i,k) if its fluctuations about the nominal values k produced sensitivity coeffi- 
cients or Kolmogorov-Smirnov coefficients larger than 10%. Equivalently the state yi was deemed not 
to be sensitive to variations of those parameters kj for which the Taylor sensitivity coefficients or the 
Kolmogorov-Smirnov sensitivity coefficients were smaller than 10%. 

The two methods are conceptually different and hence the sensitivity coefficients Sjj (ti , k) and S.y (k) 
introduced in ([3]) and (|7J) respectively, are not expected to be numerically similar to the Kolmogorov- 
Smirnov coefficients relative to the same processes. Nevertheless they exhibit similar qualitative results, 
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in the sense that most of the parameters that are important by the Taylor expansion method, are 
likewise important by the Monte Carlo filtering method. We first discarded those parameters to which 
neither Ca 2+ and RhoGTP are sensitive by the quantitative criterion indicated above. Then in Table |5] 
we cross-listed those parameters to which Ca 2+ or RhoGTP or both were sensitive. According to the 
adopted criterion of importance, our analysis shows that kj and kg are the most important parameters 
to reproduce the expected system behavior, because they respectively have ten and six entries in Table [8] 
In Figs. [2] and [3] , as an example, we show the cumulative distribution functions (fit))) calculated by Monte 
Carlo filtering for the parameter k-j in the case j/j = Ca 2+ and y t = RhoGTP respectively. The red curves 
represent the distribution of "acceptable" parameters, whereas the blue ones represent the distributions of 
"unacceptable" parameters. In panel A we report the results relative to the sensitivity test performed with 
the objective function calculated as in (fT2|) . panel B shows the results of the Monte Carlo filtering method 
with the objective function estimated according to (|13|) . and panel C and D are relative to the choice 
of the objective functions in (| 16[) respectively. Similarly, In Fig. 01 we show the cumulative distribution 
functions (jTU)) calculated by Monte Carlo filtering for the parameter &26 m the case Vi = RhoGTP. 
Examination of Tables [THE] permits one to classify the input parameters into three broad categories: 

RhoGTP-Only Sensitivity Parameters (four or more entries in Table |8|) 

lt26 PARI to G13 binding rate; 
k 36 G13GTP-GEF dissociation rate; 
k 37 GEF hydrolyzing G13GTP; 
k 40 GEFG13* to Rho binding rate; 
k 4 3 Rho deactivation rate. 
Variations of these parameter affect only RhoGTP and its functionals 



[RhoGTP] (ti,k), 



J [RhoGTP] (t, k)dt 



(18) 



[RhoGTP] max (k), 



20 



introduced in (U)-©, and (|17[) . Neither small variations (Taylor's expansion), nor large variations 
(Monte-Carlo filtering) of these parameters about their nominal values affect Ca 2+ and its function- 
al. This is an expected result as these parameters only belong to the RhoGTP module downstream of 
G13 (see Fig. [J). 

Interestingly, our sensitivity analysis found, independently of j4], that k 2 6, which represents PARI 
binding rate to G12/13, is a key factor in this model to reproduce the experimental data. This emerges 
from two different sensitivity analysis methods and constitutes a further test of validity for the signaling 
model adopted in [4]. 

Ca 2+ -Only Sensitivity Parameters (four or more entries in Table [8]) 

ki 9 rate of hydrolysis of GqGTP by PLC/3; 
k 22 PLC/?-PIP2 binding rate; 
k 2 5 IP3 rate of consumption; 
Variations of these parameter affect only Ca 2+ and its functionals 

[Ca 2 +](i 4 ,k), / T [Ca 2 +](i,k)d< 

(19) 

[Ca 2 +] max (k), t£*(k) 

introduced in (U])-©, and (fTT)) . Small variations (Taylor's expansion) or large variations (Monte-Carlo 
filtering) of these parameters about their nominal values do not affect RhoGTP or its functionals. Again, 
this is an expected result as these parameters only belong to the Ca 2+ signaling module downstream Gq. 

The presence of parameters that affect the Ca 2+ output but not the RhoGTP output and viceversa 
is due to the modular structure of the model in [J]. In this model, Gq and G13 pathways are described 
by two separate computational modules with a common input, i.e. the activated PARI (Fig. [I}. 

The Ca 2+ functionals in (|19p are not sensitive to small variations of k 2 5 (rate of IP3 consumption) 
about its nominal value k 2 5 (Taylor's method for Ak 2 5 < 5%). They are however severely sensitive for 
variations of k 2 5 in the range (10 _1 k 2 5, 10k 2 s) (Monte Carlo filtering method). 

In all cases however, the notion of "sensitivity" and "relevance" is the same by both methods, pointing 
to a robustness and self-consistency of the model in [3J. 



Table 8. List of the most important parameters as deemed by the six different sensitivity tests adopted 





Sensitivity Anlysis Method 






Taylor U 


Filtering ti 


Taylor T 


Filtering T 


Filtering y m ax 


Filtering t maT 




k 


Ca 2 + 


RhoGTP 


Ca 2 + 


RhoGTP 


Ca 2 + 


RhoGTP 


Ca 2 + 


RhoGTP 


Ca 2 + 


RhoGTP 


Ca 2 + 


RhoGTP 


FUNCTIONAL ROLE 


3 






X 












X 


X 






Consumption PAR.l-AP 


7 




X 


X 


X 




X 


X 


X 


X 


X 


X 


X 


PAR.1 deactivation 


8 








X 






X 




X 


X 


X 


X 


PARI binding Gq 


9 
























X 




14 






X 












X 




X 


X 


Gq releasing /?7 
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X 








X 




X 




X 




PLC£ hydrolyzing GqGTP 


22 






X 








X 




X 




X 




PLC0 binding PIP2 


25 






X 








X 




X 




X 




Consumption of IP3 


26 




X 




X 




X 




X 




X 






PARI binding G13 


27 








X 












X 








34 












X 












X 


G13 deactivation 


35 




X 




X 












X 






GEF binding G13GTP 


36 




X 




X 








X 




X 








37 








X 




X 




X 




X 




X 


GEF hydrolizing G13GTP 


40 




X 




X 




X 




X 




X 






GEFG13* binding Rho 


41 








X 




















43 




X 




X 




X 




X 




X 






Rho deactivation 



Comparison of the sensitivity analysis results for the most important parameters reported in bold font in Tables HHT] The symbols x 
mark the importance of a given parameter kj according to the different sensitivity test. Parameters k*i and k% are important in the 
majority of the sensitivity tests performed. 
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Ca 2+ -RhoGTP Sensitivity Parameters 

kj PARI deactivation rate; 

k 8 PARl-Gq binding rate; 

ki4 rate of release of (3j by Gq; 

PARI acts on the common part of the Ca 2+ and RhoGTP pathways and accordingly both of them and 
their functionals are sensitive to both these parameters. While PARI deactivation rate k7 is important 
to both Ca 2+ and RhoGTP, the various components of the system respond differently to variations of 
this parameter. First, Ca 2+ and its functionals (TK))) . are not sensitive to small variations of about its 
nominal value kj (Taylor's methods for \£k?\ < 5%). On the other hand all the Ca 2+ functionals are 
very sensitive for variations of k7 in the range (10 _1 k7, lOk.7) (Monte Carlo filtering method). 

Contrarily RhoGTP and its functionals are very sensitive to any variation k7, small or large, and by 
whatever method sensitivity is evaluated (Taylor or Monte Carlo filtering). 

Finally, we observe that the pathways analyzed in this work are shared in a number of signaling 
cellular systems. In a recent work, we devised a mathematical model of PARl-mediated signaling in 
human platelets [2] that shares several features with the model in [4]. For these reasons and for the 
generality of the methods adopted in this work, the conclusions of this study can, at least qualitatively, 
be extended to the model in [2]. 
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